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We investigate the linear propagation of a paraxial optical beam in anisotropic 
media. We start from the eigenmode solution of the plane wave in the media, 
then subsequently derive the wave equation for the beam propagating along 
a general direction except the optic axes. The wave equation has a term 
containing the second mixed partial derivative which originates from the 
anisotropy, and this term can result in the rotation of the beam spot. The 
rotation effect is investigated by solving analytically the wave equation with 
an initial elliptical Gaussian beam for both uniaxial and biaxial media. 
For both media, it is found that there exists a specific direction, which is 
dependent on anisotropy of the media, on the cross-section perpendicular 
to propagation direction to determine the rotation of the beam spot. When 
the major axis of the elliptical spot of the input beam is parallel to or 
perpendicular to the specific direction, the beam spot will not rotate during 
propagation, otherwise, it will rotate with the direction and the velocity 
determined by input parameters of the beam, (c) 2010 Optical Society of 
America 

OCIS codes: 260.1180, 260.1960. 
1. Introduction 

The propagation of light in anisotropic media is still receiving a good deal of attention 
from both the experimental and the theoretical points of view. The subject is of funda- 
mental interest because of its important role in polarizers, compensators, amplitude- and 
phase- modulation devices, also, in nonlinear optical phenomena which usually take place in 
anisotropic media. An accurate electromagnetic description is required for the evaluation of 
the field distribution in, for example, designing polarizers and compensators or in describing 
the nonlinear processes like second harmonic generation. 
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A number of studies [1-12] have been devoted to an analysis of optical beam propagation in 
anisotropic medium, but most of them are concerned about the uniaxial cases. For instance, 
Fleck et al. [1] derived the paraxial equations for both the ordinary and the extraordinary 
components of the field in uniform uniaxial anisotropic media. The derivation is based on 
the paraxial approximation that the electric field corresponding to the beam is plane po- 
larized with a polarization corresponding to either an ordinary or an extraordinary wave. 
Ciattoni et al. [2-6] proposed a theoretical scheme that the exact plane wave angular spec- 
trum representation is approximated within the paraxial constraint, obtaining the paraxial 
expressions for the ordinary and extraordinary components of a field propagating in uniaxial 
media. Trippenbach et al. [7, 8] considered the propagation of light pulses in anisotropic and 
dispersive media, deriving the equation for the slowly varying envelope of the electric field. 
The equation has applications in divers fields. In their paper, it was used to investigate the 
propagation of an optical pulse in uniaxial media, and found that the equation contains 
terms that rotate the pulse in time domain. The rotation effect that originate because of the 
anisotropy of the medium and the finite size of the field. The similar effect in spatial domain 
can be found in the beam propagation in anisotropic media, however, it remains undiscussed. 

Generally speaking, investigation of optical beam propagation in biaxial media is more 
complicated than that in uniaxial media because of the lower symmetry in the dielectric 
tensor of the former. To our knowledge, the first attempt to quantify the propagation of 
beam in biaxial media was made by Dreger [9]. In his paper, he derived the optical beam 
propagators for biaxial media and specialized them for three cases of wave- vector direction: 
general directions, coincident with an optic axis, and nearly parallel to an optic axis. His 
derivation refers to complicated mathematical analysis. In our paper, we will give another 
way of deriving the wave equation for paraxial beam propagating in the general directions. 

The paper is organized as follow: In section 2, we use the method that has been introduced 
in [13] to obtain the eigenmode solution in biaxial media. This method is well known, but to 
our knowledge, the results for biaxial media are unavailable in the literatures published so far. 
In section 3 we derive the equation for the optical beam in anisotropic media. Previously, Qi 
Guo and Sien Chi [10] have derived the nonlinear paraxial wave equation for the propagation 
of an optical beam in anisotropic media. Because the nonlinear effect are treated as a small 
perturbation, the whole derivation is available for the linear situation if the nonlinear term 
is not considered. Besides, we find that there is an improper approximation in the process 
of the derivation. In this section, therefore, we adopt their way to derive the beam equation 
and make an improvement in the derivation. In section 4 we illustrate the dynamics of the 
beam propagation by using initial elliptical Gaussian beams. Section 5 is a summary. 
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2. Eigenmodes in anisotropic media 

We start from the eigenanalysis of optical field in anisotropic media. From the beginning to 
Eq. (14) in this section, we briefly recall the principles of the methods of how to obtain the 
eigenmodes, referring the readers to previous publications for details. 

For a monochromatic optical field E(r, t) = E(r) exp(zci;t) in a magnetically uniform 
anisotropic medium, E(r) is described by the wave equation [1] 



V X [V X E(r)] - "^e ■ E(r) = 0, (1) 



where e is the relative dielectric tensor. Assuming that Eq. (1) has a plane wave solution 

E(r) = Eoexp(-zk-r), (2) 

where the amplitude Eq is a constant and k is the wave vector, we have the following 
eigenvalue equation [14] 

(kk - /c^J + —e) ■ Eo = 0, (3) 
where / is an unit tensor of rank 2. Moreover, the dispersion relation 

dei{\^\i-k'^i +—e) = Q (4) 

can be obtained from the necessary condition for nontrivial solutions to Eq. (3). On the other 
hand, the electric field E and the electric displacement D are related by the constitutive law, 

D = ■ E, (5) 

where Eq is the dielectric constant in vacuum. The Gauss's law for the plane wave requires 
that 

k-D = 0, (6) 

implying that k is perpendicular to D. It is proved [13, 15] that for every direction of prop- 
agation k/|k| for the plane wave in anisotropic media, there are, in general, two orthogo- 
nal polarizations of D and, correspondingly, two values of k. The allowed plane waves are 
called eigenmodes and their corresponding polarization directions D/|D| are called eigen- 
vectors [15]. As soon as D/|D| are found, the polarizations of E can be obtained from the 
inverse function of Eq. (5), 

D 

(7) 



E _ - ^ 



lEi i£-i-^r 

1^1 1^ |D|I 

where is the inverse of the tenser e. One convenient way of finding D/|D| can be visualized 
with the aid of the index ellipsoid [13], so we are going to solve the eigenmodes problem in 
biaxial media by this way. 
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Consider a lossless anisotropic medium, in the principle coordinate system (x', y', z'), char- 
acterized by primes, the dielectric tensor e' is of a diagonal form, 

/ Er' \ 



e' 



ey> 



(8) 

V e^, j 

where Ex'i Syi and s-^ are real, and e^i > £y' > £z' is supposed by convenience. The equation 
of the index ellipsoid is [13] 

'-+y-+'- = i, (9) 

£x' £y' £z' 

as shown in Fig. 1(a). We write the wave vector k of the plane wave in the form of 

k = Gx'/cx' + Gy'ky' + e^'kz', (10) 

where e^r, ey/ and e^,' are unit vectors along x', y' and z' axes, respectively, while kx' = k cos a, 
kyi = k cos P and kz' = k cos 7, with cos a, cos (3 and cos 7 the direction cosines. The equation 
of the plane S perpendicular to k through the origin is 

a;'cosa + ?/'cos/3 + z'cos7 = 0. (11) 

This plane cuts the ellipsoid in the ellipse (Fig. 1 (b)) from which the eigenvectors and the 
values of k can be found [13]: 1, the magnitudes of the major and minor axes represent the 
refractive indices of the eigenmodes which are related to k by /c = un/c, where c is the phase 
velocity in vacuum; 2, the eigenvectors are along the corresponding major and minor axes, 
respectively. 

Since the magnitudes of the major and minor axes, denoted as ni and 77,2 as shown in 
Fig. 1(b), are the longest and shortest diameters of the ellipse, they can be obtained by the 
Lagrange multiplier method. We can determine rii and n2 by finding the extremum of [13] 

subject to the conditions Eq. (9) and Eq. (11). According to the Lagrange multiplier method, 
we construct the function [13] 

f = + y'^ + z'^ + XJ — + ^ + — - 

ySx' Eyi Sz' 

+A2(x'cosa + y'cos/3 + 2'cos7), (13) 

where Ai and A2 are the two undetermined multipliers. Now the problem is equivalent to 
finding the extremum of / subject to no subsidiary conditions. The necessary conditions for 
the extremum of / are that its derivatives with respect to x, y and z should vanish [13], i.e. 

-^ = 2x' + i + A2Cosa = 14a 

ax Syi 
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— = 2y A h A2 cos = 

oy Syi 

df , 2z'\i , 

-^ = 2z' + + A2 COS7 = 0. 

OZ 6,1 



(14b) 



(14c) 



By solving the three equations (9), (11) and (14), we can obtain the solutions of x', ?/', z' , 
Ai and A2, then the two extremum of n (n only takes the positive values) are 



2a 



(15) 



where j = 1,2 and 



a 
b 



cos^ a + eyi cos^ /3 + e-^i cos^ 7 
(e^,/ + )e^/ cos^ 7 + (e^/ + e^/)ey/ cos^ /3 + (ey + e^/)e3;/ cos^ a. 



The plus sign is associated with ni, the minus sign is associated with n2, respectively. The 
values of k are given by kj = rijU / c. We call the eigenmode propagating with the wave vector 
ki the mode 1 wave, and that with k2 the mode 2 wave. With the solutions of Ai and A2 
inserted into Eq. (14), we get the eigenvectors: 



ID. 



e^' cos a eyi cos /3 e^' cos 7 



By 



ey. - n- 



1 I ' 



(16) 



where 



ex' cos a 
er' — 



+ 



eyi COS (3 



It is easy to prove that Di ■ D2 = which means that the two eigenvectors are mutually 
orthogonal to each other. 

We can see that rti = 77-2 if the polynomial under the radical sign in Eq. (15) is equal to 
zero, 



ez' cos 7 



Aaex'ey/ez' = 0. 



(17) 



As Ex' > Syl > Ez', it can be rewritten as 



[ex'iey' - ez>)cos a - (e^./ - e^^/) cos 7] + ey,{ex' - e^') cos /3 + cos /3 
x{2ey,{ex' - ez')[ez'{ex' - e^/) cos^7 + ex'{ey' - e^') cos^a]} = 0, 

where the relation that cos^ a + cos^ /3 + cos^ 7 = 1 is used, then we can obtain 

cos /3 



cos^ 7 



0, 



cos^a ez'iex' - eyi) 



tan ar- 



(19a) 
(19b) 
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Two directions are specified implicitly by Eq. (19). Along these two directions, we will have 
rti = 722 which shows that the ellipse in Fig. 1(b) is a circle, then the whole formulation 
breaks down since the Lagrange multiplier method is unavailable in such case. These two 
special directions are known as the optic axes [13]. Eq. (19) shows that they lie in the x' — z' 
plane, at angles ±ac to the x'-axis. As the optic axes are the singularities of the dispersion 
relations [1], in the next section about beam propagation, we only discuss the cases that the 
optical beam is neither directed along an optic axis nor so close to an optic axis that any 
(spatial) spectral component of the beam is directed along an optic axis. The results for the 
uniaxial case can be obtain by assuming that e^' = ^y' or eyi = e^/, since these results are 
available in many literatures [13], we don't present them here. 

3. Wave equation for paraxial beam in anisotropic media 

The wave equation for paraxial beams in anisotropic media had been derived in reference [10], 
but the derivation has an improper approximation in the process. In the first subsection, we 
will follow the same approach to derive the wave equation and make an improvement. The 
wave equation contains coefficients that account for different propagation behaviors. The 
expressions of these coefficients are needed to quantify the beam propagation, and they have 
been obtained in the case of uniaxial media [10]. In the second subsection, we will discuss 
the case in biaxial media and find the expressions of these coefficients. 

3. A. Derivation of wave equation for paraxial beams 

Now, we consider the propagation of an optical beam, which is formulated by a group of 
plane waves that belong to only one of the two mutually orthogonal eigenmodes. The solution 
of 1, therefore, can be expressed as 



where k = K + ez/c^, r = R + e^^ (ez is an unit vector along 2-axis, the capital letters K and 
R represent transverse wave vector and transverse coordinate vector perpendicular to e^, 
respectively, K = e^k^ + Gyky, R = e^x + eyy, with ex and ey defined similar to e^). In the 
paraxial beam C3iSG, clS the wavevectors of the component waves do not deviate much from 
the central wave vector kg, it is convenient to introduce a new coordinate system {x,y,z), 
termed propagation coordinate system hereinafter, so that its 2;-coordinate axis coincides 
with kg. In this coordinate system, the input field in z = reads 



In the spirit of the paraxial approximation, we can assume that the electric field corre- 
sponding to the beam is linearly polarized with a polarization corresponding to the eigenmode 




(20) 




(21) 
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propagating with the center wave- vector kg [1]. Concretely, if we write the amphtude of the 
component waves as Eo(K) = p(K)£'o(K), where p(K) is a unit vector representing the 
polarization state, according to the approximation, we have 

p(K) ^ p(0) = po. (22) 

This approximation is not the same as that shown in [10], where the whole amplitude Eo(K) 
is approximated as Eo(K) ^ po-Eo(O). This approximation is improper because Eo{0) is a 
constant as soon as the input field is specified, thus it will make Eg. (21) unable to describe 
an arbitrary input field. 

under the paraxial condition, kz can be expanded in a Taylor series at the point K = 0, 

h{K) = ko + VKko ■ K + ^WK^Kko : KK + • • • , (23) 

with Vxko = ^Kkz\K=o and so on, and V^- = d/dk^e^ + d/dkyCy. If we keep terms up to 
second order in Eg. (23) and make use of Eg. (22), we can obtain from Eg. (20), 

E{r) = poA{r)exp{ikoz), (24) 

where 



A{r)= / / Eo(K)exp[2<l>(K,R,z)]dK (25) 

J J — oo 

is the slowly varying amplitude, and the phase factor $ reads 

$(K, R, z) = (Va'A;o ■ K + ■ KK)z + K ■ R. 

Introducing Eg. (24) into Eg. (1) with the help of Eg. (3) and neglecting the second derivative 
d'^A/dz^, we can obtain [10] 

{kz - A;o)(Ke^ + e^K) • Eoexp(z<l>)dK + (VrG^ + g^Vr) ■ po^ 

oz 




(e^e^ - /) ■ Po 



dA 

2tko-^ + / / Eoik'i - e,) exp(z$)dK 




0. (26) 



with = d/dxGx + d/dyGy. It was proved that the first two terms of Eg. (26) can be 
neglected [10]. From the expansion (23) we make the approximation for the factor [k1 — k^) 
in the integration, 

kl-kl^ 2ko{kz - ko) = 2koVKko ■ K + k^V xk^ : KK. (27) 

Substituting Eg. (27) into Eg. (26), we obtain the beam propagation equation [10] 
.,dA , dA ^ dA^ 1 / d^A d^A ^ d^A\ ^ ^ 
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where 5i = dk^ko, 



X or y. Eg. (28) is a general form of the evolution 



equation for the beam in anisotropic media. 6i and 6ij are named 6 coefficients hereafter, 
and each one has its physical meaning. The coefficients 6^^ and 6y describe the walk-off of the 
beam center in x and y directions, respectively. The coefficients S^x and 6yy are the general 
Fresnel diffraction coefficients which are modified by the anisotropy of the medium and are 
in general not equal. The coefficient S^y involving the mixed derivative d'^A/dxdy vanishes in 
isotropic media; it is nonvanishing only if the index of the refraction depends on the direction 
of propagation; it can enhance the diffraction and rotate the cross-section of the beam in the 
X — y plane. These effects will be illustrated later. 

It should be reminded that the field is determined by Eg. (24) in the lowest-order approxi- 
mation. Since E(r) is proportional to poA(r), the vectorial field is completely specified when 
the scalar beam equation Eg. (28) is solved for A{r). To find the concrete expressions for dif- 
ferent modes, we should first obtain po and its accompanying function kz(K) by making use 
of the equations Eg. (16), Eg. (4) and Eg. (15), and then obtain ko and the 6 coefficients. In 
the next section, we will discuss how to obtain the 6 coefficients in a general biaxial medium, 
and how to reduce them in the uniaxial medium case. 



3.B. Expressions of 6 coefficients for biaxial media 

Generally, available datas of the dielectric tensor are given in the principle coordinate system 
{x',y',z') [16], in which i' is of a diagonal form, but the paraxial equation Eg. (28) for 
the optical beam is obtained in the propagation coordinate system {x,y,z). To discuss the 
problem, therefore, the dielectric tensor should be transformed from the principal coordinate 
system to the propagation coordinate system, in which the dielectric tensor is denoted by e. 
To do this, we use the so-called " x-convention" definition of the Euler angles {(f),6,ip) [17], 
as shown in Fig. 2. In this convention, the coordinate rotation can be written as 



(29) 





X 


\ 




/ 


x' 


\ 




y 




= T 




y' 




v 


z 


) 




V 


z' 


I 



with the rotation matrix T, 



/ cos il) cos — cos Q sin sin -0 cos -0 sin -|- cos d cos sin -0 sin -0 sin Q 



\ 



sin 0^ cos — cos Q sin cos cos Q cos cos ^/^ — sin -0 sin cos sin 6 
sin Q sin — sin Q cos cos Q 



(30) 



Through the angles and 6*, we can rotate the z' -coordinate axis to the z -coordinate axis, 
i.e. the direction of the central wave vector ko. In the principle coordinate system, a given 
propagation direction ko/|ko| is described by the direction cosines (cos a, cos /3, cos 7), then 
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the values of Euler angles and 6 can be obtained from the relation 



e = J (31a) 

cos a , 

tan0 = -. (-jlb) 

cos p 

The nutation angle tp of the final rotation can be arbitrary chosen since it is used for the 
orientation of the x and y axes. In general, it is chosen to align the x and y coordinates with 
the two orthogonal eigenvectors, and these special values of ip, satisfy 

-(a cot + tan 0) 

tan^o = 7z ^ ^ — , 32 

(1 — a) cosy 

with a = {nl — ex')£y'[{nl — eyf)ex']~^- The angle ipo is to be distinguished from ip which 
can be chosen arbitrarily in defining the propagation coordinate system. If ip = ipo, we have 
Di/|Di| along x-axis and D2/ID2I along y-axis, respectively. As a result, ipQ can be used to 
denote the directions of the two eigenvectors in propagation coordinate system. 

With the coordinates rotation, the dielectric tensor in the propagation coordinate system 
becomes 



^11 ^12 ^13 

i = Te'T"^ = I 622 £23 



\ 



(33) 



£^13 £^23 £33 I 

where T"-*- is the inverse matrix of T, and the components of e are, respectively. 



£11 = (cos'^/'sin^ — cos^sin'?/'sin0)^e2;, 

+ (cos Q cos sin + cos -^p sin (jyfeyi + sin^ 6 sin^ ipEz' , 
£12 = (— cos 0sin'?/' — cos6' COS'?/' sin0)(cos'?/; cos — cos^ sin-?/' sin 0)63;/ 
+ (cos 9 cos sin ip + cos ip sin 0) (cos 9 cos ip cos — sin ip sin (p)eyi 
+ cos Ip sin^ 9 sin ipEzi , 
£13 = (cos '0 cos — cos ^ sin -0 sin 0) sin ^ sin 0ea;/ 

— cos sin ^(cos 9 cos sin ip + cos ip sin (p)eyi + cos 9 sin 9 sin ipSz' 
^22 = (— cos sin?/; — cos 6* cos 0^ sin 0)^62,/ 

+ (cos 9 cos Ip cos — sin -0 sin (pYsyi + cos^ ip sin^ ^£2', 
^23 = (cos sin -0 — cos 9 cos ip sin 0) sin 9 sin 0e2;/ 

— cos sin ^(cos 9 cos ?/; cos — sin ip sin 0)5^' + cos 9 cos ?/; sin 9ez' 



£33 = sin^ 9 siv? (pEx' + cos^ sin^ + cos^ 9ez'. 



Substituting Eg. (33) into Eg. (4), we get the dispersion relationship in propagation coor- 
dinate system. The expressions of 5 coefficients can be obtained from the procedure that 
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introduced in the last paragraph in section 3. A. Since the calculation is complicated, we just 
present the results here: 

,2 



^13^22 — ^12^23 — ^13'^,- 



'^e-SSn^i - ^11^33 - ^22^33 + + g 



23 



^11^23 ~ ^12^13 — g23'^i 



2£:33'^j - ^11^33 



^22^33 + g 



13 



-2 ' 
-23 



1 

ko 



ko 



-^12 - ^13 + ^11 (^22 + £33) - ny^n + £33) - 4:n^j5xi£i3 + £33Sx) 



2g33'T'f - ^11^33 - ^22^33 + sj^ + E 



+ 61 



23 



-12 



-23 



^22(^11 + £33) - n^j[e22 + ^33 



Anpy{e 23 



^336y) 



2es3nj 



Jxy 



1 

ko 



^12^33 ~ ^13^23 



n]ei2 



^11^33 "~ ^22^33 + ^13 + ^23 
2^ X 0^2, 



+ 61 



2nje23Sx - 2njei3Sy - Anje33SxS^ 



2^33^^^ ~ ^11^33 ~ ^22^33 + ^13 + ^23 



+ SxSy 



(34a) 

(34b) 
(34c) 
(34d) 
(34e) 



where ko = riju/c with rij given by Eg. (15). If Ui is substituted, we will obtain a set of 
6 coefficients for the mode 1 beam and that for the mode 2 beam with n2 substituted. In 
particular, ip = ipo can reduce the 6 coefficients to the same forms that presented in [9]. 

Eg. (28) together with the 6 coefficients describes the paraxial beam propagating in a 
general direction except along the optics axes in a general anisotropic medium. In the uniaxial 
case, for example, the positive uniaxial crystal with Ex' = Syi = and e^/ = n"^ in e', the 
mode 1 beam turns out to be the ordinary beam with ni reduced to Uo, and the S coefficients 



are reduced to 6r 



yy 



-c/{ujno) and zero of the others, then we obtain the standard 



paraxial equation for the ordinary beam. 



.dA c 
dz 2ojnn 



d^A d^A 



A = 0. 



(35) 



On the other hand, the mode 2 beam is the extraordinary beam with the 5 coefficients 
reduced to 

{ni — ni) sin ip sin 29 



^ 2{n'l siv? + nl cos^ Oy 

{n1 — n1) cos 'ip sm.26 
^ ~ 2{nl sin^d + nl cos^^)' 

1 ra^ cos^ tp sir? 6 + ninl (sin^ ip sin^ 6 + cos 



k 



ini sin 6 



n1 cos^ I 



1 rig sin ipsm. 6 + ninlicos ip sm. 6 + cos 



^yy 



Jxy 



{n1 sin^ 9 + n1 cos^ 9Y 
1 n^(n^ — Wg) sin2'?/'sin^ ^ 

k^^) 2{nlsiYi^ e + ni cos^ey ' 



(36a) 
(36b) 
(36c) 

(36d) 
(36e) 
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where k^^ = n^'^'^uj/c and 



is the refractive index with the form expressed in the propagation coordinate system. Eg. 
(37) is obtained by making the coordinate transformation on n2 after it has been reduced 
under the uniaxial condition. Here we pay some attention on the results for the extraordinary 
beam. The symmetry about the 2;'-axis of the index ellipsoid makes the propagation direction 
only described by ^, so the 5 coefficients are independent on the angle 0. For the coefficient 
5xyi Eg. (36e) shows that once 6 is specified, its value is modified by the angle ip and the 
sine function makes the value periodic. To physically understand these properties, we recall 
that the wave vectors of the extraordinary plane waves belong to an extraordinary eUipsoid 
of semi axes UoUj/c, n(.u/c and n(.uj/c, of which the surface is called the normal surface [15]. 
For the extraordinary beam, the curved surface of the ellipsoid near the central wave vector 
kp*^^ is not rotational symmetry about kg'^'' if kg^'* is not along the optical axis. Therefore, 
the mixed derivatives terms arise in the Taylor expansion Eg. (23) for a general ip, and its 
periodic characteristic is readily understood by recalling that the orientation of the x and 
y axes repeat again if ip increase from zero to 2tx. The coefficient 6xy in the case of biaxial 
media has the similar properties but it is not evident from Eg. (34e). 

In this uniaxial case. Eg. (19) implies that the optical axis coincides with the 2;'-axis. If we 
set = 7r/2 which lies the x-axis in the plane that formed by the optic axis and the 2;- axis, 
and substitute (38) into Eg. (28), we will get the wave equation for the extraordinary beam 
which is equal to the linear part of the nonlinear equation obtained by Qi Guo and Sien Chi 
(Eq. (42) in [10]), 



dA (n2-n2)sin2e dA 

dz 2{n'lsm^ 9 + nl cos^ 9) dx 
1 d^A^ 



nl d^A 



{nl sin 9 + nl cos^ 9)"^ dx 



nl sin 9 + nl cos^ 9 dy^ 



0. (3^ 



Further more, 9 = 7t/2 will reduce Eg. (38) to the form of describing the paraxial beam 
propagation orthogonal to the optical axis, which is the same result given by the equation 
(8) in reference [6]. 

4. Elliptical Gaussian beam propagation 

The properties of paraxial beam propagation in anisotropic media will be fully elucidated by 
using a specified input field. In the first subsection, we will use an initial elliptic Gaussian 
field for investigation and obtain the analytic solution. We find that the elliptic Gaussian 
beam will rotate during propagation and the rotation is discussed in detail. In the second 
subsection, we will give some examples for illustration. 
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4- A. Analytic solution of initial elliptic Gaussian beam 

We discuss the elliptic Gaussian beam propagation in anisotropic media. Suppose that the 
input field at 2 = is given by 

Eo(x, y, 0) = poAo exp (^^^ exp (^^^ , (39) 

where Wx and Wy are the initial beam widths in x and y directions, respectively, and for 
convenience, Wx > Wy] Aq is the amplitude of the field. With this input field, the analytic 
solution of Eg. (28) is [7] 



exD < -\-^^^^-y^^lH-^^yy^y^- . 



_y2 



^[Wl - iSxxZ + 6lyZ^/{wl - i6yyZ)]{wl - l6yyZ) 

where X = x + SxZ, Y = y + 6yZ, respectively. To investigate the field evolution, we should 
first obtain the modulus of the complex expression Eg. (40), which reads 

M exp[-(§X2 + §r2 + fXF)] 
\A{x,y,z)\ = ^ ^ ^, (41) 

where 

Bi = W^Wy + Z Wjyy + Z Wyd^y, 

D 42|22c-2|22c-2 

B2 = W^Wy + Z Wyd^^ + Z Wj^y, 

B3 = -z^6xy{wl6xx + wl6yy), 

B = Z\wl6yy + Wl6xxf + \wlwl + Z^{bly - bxx^yy)f- 

The quadratic form of the exponent shows that the field distribution is a rotated ellipse in 
X — y plane if -B3 7^ 0. i?3 is proportional to 6xy^ so 8xy = for no rotation. The rotation 
angle (p, defined as the tilt of the major axis relative to X-axis, satisfies 

B3 



tan2(/? — 

Bi — B2 

2z'^Sxy{wlSxX + wlSyy) ^^^^ 



A word of conventions is needed here. The rotation angle ip determined by Eg. (42) has a 
range within ±7r/2, so it will be the convention of this paper to fix the quantity unambigu- 
ously by defining the arctangent function as having a range lying in the interval [— 7r,+7r]. 
On this convention, the major to semimajor axes ratio 77 of the rotated ellipse is 

Sin^ (y9 + i?2 COS^ (/3 — ii?3 sin 2<y9\ ^''^ 

^ ~ V^i cos2 P + B2 sin^ if + \B^ sin 2p ) ' ^ ' 
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When 5xy 7^ 0, the rotation velocity, described by dip/dz, can be obtained from the deriva- 
tive of Eg. (42) with respect to z, 



dip ^ 2cos\2ip)zwlw^{wl - w^)5xy{wl5xx + w^Syy) 



(44) 



Generally, the diffraction coefficients d^x and 5yy are negatives, making the beam spread in 
the medium. As a result, the rotation direction is determined by the coefficient 5xy which is 
directly associated with the angle ip while the propagation direction has been specified. For a 
chosen if), if 5xy < 0, we have dip/dz > 0, implying that the beam will rotate counterclockwise. 
On the contrary, 6xy > is corresponding to the clockwise rotation of the beam. For large 
z, Eg. (42) is approximate to 



2A + 



yy 



tan2(/3 ^ 5 1 . (45) 

A2 _ !ii^2 I |'!£i_1^A2 

«xx u,l"yy^ywl ^Pxy 

Eg. (45) shows that the rotation angle will approach a value that depends on the ratio of 
the initial beam widths Wx/wy. As a special case, when the input field is a circular Gaussian 
beam, i.e. Wx = Wy, it can be proved directly from Eg. (40) that the input field will lose its 
circular symmetry during propagation and become elliptical. This is because the physical 
properties of the diffractive spreading (described by 6xx, ^yy, ^xy) are different in every 
transverse directions in an anisotropic media. With Wx = Wy, Eg. (42) gives 

tan2(^ = -^^. (46) 

(^xx "yy 

The rotation angle is independent on z, so the beam will not rotate but incline by an unvaried 
angle during propagation. 

For the case of 5xy = 0, Eg. (40) can be factorized into the product of two one- dimensional 
Gaussian beams. In such case, the beam will not rotate but just spread in x and y directions 
with the diffraction coefficients 5xx and 5yy, respectively. 

Here we should remind that the angle ip can be freely chosen and plays a role of fixing the 
X and y axes. Since the input field Eg. (39) is of a regular form in the propagation coordinate 
system, different values of ip correspond to different input fields. If we choose a value of ip 
that makes 5xy 7^ 0, then the input elliptic Gaussian beam will rotate during propagation, 
otherwise, there will be no rotation effect on the beam. To illuminate the relationship between 
ip and the input field, we still take the extraordinary beam in the positive uniaxial crystal 
for example. The 5 coefficients for the beam are given by Eq. (36). We can obtain from Eg. 
(36e) that 5xy = Q ii ip = m7r/2, where m = 0, 1, 2, 3. If ip takes one of these values, which 
corresponds to lying one of the waists {wx or Wy) in the plane that formed by the optic axis 
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and ko, see Fig. 3, the extraordinary beam will propagate with no rotation effect (as the 
case of = 7r/2). Otherwise, the extraordinary beam will rotate during propagation (as the 
case of = vr/4). The cases of the biaxial media are similar except that the values of for 
6xy = are different from that in the uniaxial cases. 

4.B. Examples 

Here we will use some examples to illustrate the properties of the propagation that discussed 
above. For more general, we show the beams propagation in the biaxial crystal KNbOs. The 
wavelength of the input beams is A = 514nm, and for this wavelength, the elements of the 
dielectric tensor in principle coordinate system are reported as e^' = 5.736, Ey' = 5.446, 
Ez' = 4.893 [16] which place the optic axes at ac = ±56.21°. The optical beam is supposed 
to propagate along a general direction (a, f3, 7) = (84.36°, 88°, 6°) in the principle coordinate 
system. This direction is arbitrary chosen but should not be too close to the optic axes. It 
will be convenient to use this direction in all examples. The refractive indices of the two 
eigenmodes are rii = 2.393 and 77-2 = 2.334. According to Eg. (31), the first two Euler angles 
are (f) = 109.50°, 6 = 6.00°, and from Eg. (32), we have tpo = 69.94° . Fig. 4 shows the 
parameter 6xy versus the angle ip in the supposed propagation direction. As shown, there 
are four roots for 6xy = which are 69.48° + m90°, with m = 0, 1,2,3. To accentuate the 
physics of the S^y term, we take ip = 47° in the first example to show the propagation of 
beam with S^y 7^ 0, and then ip = 69.48° is taken in the second example for comparison. 
The 6 coefficients are evaluated by (34a)- (34e), and the propagation is computed by using 
Eg. (40) with the propagation distance normalized by = zc[uj{wl + Wy)]~^. The results 
are shown in Fig. 5 and Fig. 6 in which the coordinate frames move with the centers of the 
beams, so the walk-off of the beams are not shown in the figures. 

Fig. 5 shows the propagation of the two beams of different modes with three different 
initial beam widths. Fig. 5(a) is the case of the initial circular Gaussian fields. We can 
see that after the propagation, the initially circular cross-sections of the two modes both 
become elliptical and inclined. The unvaried inclined angles given by Eg. (46) are the same, 
both are equal to 22.48°. Figures 5(b) and (c) show the propagation of the initial elliptical 
Gaussian beams. The mode 1 beams rotate counterclockwise and the mode 2 beams rotate 
clockwise, as predicted by Eg. (44). For ^ = 10, the propagation distant is long enough that 
the rotation angles can be approximately obtained by Eg. (45): the rotation angles for the 
mode 1 beams are if ~ 76.09° in (b) and ip ~ —83.65° in (c), and for the mode 2 beams, 
they are ip ^ —84.52° in (b) and if ^ —86.69° in (c), respectively. Fig. 6 shows the cases of 
Sxy = 0. As expected, there is no rotation or inclination effect on the beams. The beams just 
spread in the X and Y directions. 

To get a full picture of the beam rotation, we present the details in the following examples. 
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The parameters are the same as those in Fig. 5 except with Wx = 12.5yU.m and Wy = 10.0/im. 
The results are shown in Fig. 7 which contains four plots with different C, for each mode. 
The mode 1 beam rotates counterclockwise. The rotation angle increases with increasing 
and the cross-section keeps elliptical during propagation. The mode 2 beam has a clockwise 
rotation. For C, = 0.734, the mode 2 beam has broadened in Y due to Fresnel diffraction 6yy 
and is roughly as broad in X as in Y. Therefore, the tilt of the beam is difficult to see in 
Fig. 7(b). As the propagation proceeds, the rotation becomes evident again, as shown in Fig. 
7(c) and (d). 

Fig. 8(a) shows the rotation angles versus ^. It indicates that the rotation angles increase 
(the mode 1 beam) or decrease (the mode 2 beam) monotonously with the increasing ^, and 
for large ^, they both approach certain values which are given by Eg. (46), respectively. For 
the mode 1 beam, ip ^ 66.73°, and ip ~ —82.69° for the mode 2 beam. Fig. 8(b) shows the 
major to semimajor axes ratio rj of the ellipse describing the contour lines of \A{x,y, z)\. 
Both of the ratios share the same characteristics that each of them initially decreases due 
to the diffraction and goes through a minimum. Fig. 8(a) together with (b) shows that 
before the S^y term makes an obvious rotation on the cross-section of the beam (for example, 
ip f» 45°), rj has dramatically decreased and approached the minimum value, i.e. the ellipticity 
decreases. This is because the diffraction will exercise the most influence on the evolution of 
the beam rather than the rotation effect if the coefficient 6xy is much smaller than 6xx (or 
Syy). Although the rotation will become evident with ^ increasing, the initial rotation effect 
is strongly influenced by diffraction and is barely visible (like the mode 2 beam in Fig. 7(a) 
and (b)). As a result, in our examples we choose the spacial widths of the initial field {wx 
and Wy) to be close to each other so that the diffraction effect is not so strong to influence 
the observation of the rotation. The process of beam rotation shown in Fig. 7 may be not so 
obvious, but we think the examples are sufficient to elucidate the physics. 

An obvious beam rotation can be observed if coefficient Sxy is comparable to S^x or 6yy 
in magnitude. In order to find a large 6xy, we examine the ratio 6xy/Sxx {Syy is on the order 
of 6xx)- The ratio is a function of the the Euler angles with the maximum determined by 
{^x — ^z) /['^{^x^zY^'^] which is associated with the degree of anisotropy. According to this 
value, we find that in many natural biaxial crystals like KTP, LBO or LiNbOs [16], 6xy is at 
least one order less than 6xx- To our knowledge, natural material with large 6xy is not readily 
available yet. However, in recent scientific progress, metamaterial technology revolutionized 
modern optics and photonics by creating nearly unlimited opportunities for engineering 
material parameters. The electric and magnetic properties can be carefully designed and 
tuned by changing the geometry, size and other characteristics of meta-atoms [18,19]. As a 
result, material with large coefficient 6xy might be realized. For example, supposed that a 
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kind of metamaterial with a dielectric tenser 
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has been engineered, according to our theory above, in the supposed propagation direction, 
we have 6^^ = -8.767 x 10~^ 6yy = -4.107 x 10"^ and 6,,y = -2.339 x 10^^ for the mode 1 
beam {ip is chosen to be equal to47°). Since S^y is comparable to S^xi^yy), we can expect a 
significant effect from 6xy We present the propagation of the mode 1 beam in such material 
in Fig. 9. As expected, the inclination is obvious. 

Through these examples, we can see that the coefficient 6xy will play an rotation effect on 
the beam propagation if it is nonvanishing. Although here we just present the propagation 
of elliptical Gaussian beam, many other simulations can be made with the beam equation 
in biaxial media or uniaxial media, for example, the propagation of higher-order modes of 
Hermit-Gaussian beams, we don't present the results here. 

5. Conclusion 

We investigate the propagation of paraxial beam in anisotropic media. The propagation 
equation contains terms that account for the beam walk-off, Fresnel diffraction and beam 
rotation which is mainly discussed in this paper. The term responsible for the rotation 
originates from the finite beam size and the refractive index depending on the propagation 
direction. As a result, this term vanishes for the wave in isotropic media and for the ordinary 
wave in uniaxial media; it is nonvanishing for the extraordinary wave in uniaxial media and 
the waves in biaxial media. We investigate this rotation effect by using an initial elliptical 
Gaussian beam in the biaxial media KNbOs. The analytical solution shows that the initial 
elliptical cross-section of the beam keeps elliptical and rotates clockwise or counterclockwise 
during propagation, and for large propagation distance, the rotation angle will asymptotically 
approach a certain value which depends on the ratio of the initial beam widths. It also shows 
that an initial circular Gaussian beam will lose its circular symmetry and become elliptical 
and inclined. This is because the anisotropy makes the properties of diffractive spreading 
different in every transverse directions. Finally, we point out that the intensity of the rotation 
is associated with the degree of anisotropy. The rotation will be remarkable with high degree 
of anisotropy of the media, otherwise, it will be week and may be covered by diffraction. 
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List of Figure Captions 

Fig. 1. (a) Index ellipsoid and (b) the ellipse that the plane S cuts in the index ellipsoid. 

Fig. 2. The "x-convention" of Euler angles. The first rotation is by angle cf) about the 
z'-axis, the second rotation is by an angle 9 G [0,7r] about the former a:-axis (now oa), and 
the third rotation is by an angle if) about the z-axis. a, (3 and 7 are the direction angles of 
2;-axis in the principle coordinate system {x',y', z'). 

Fig. 3. Different values of correspond to different input fields for the extraordinary 
beam in uniaxial crystal. For ip = ir/A, the input field Eg. (39) corresponds to the dashed 
ellipse and this beam will rotate during propagation; for ip = n/2, the input field Eg. (39) 
corresponds to the solid ellipse and this beam will not rotate during propagation. 

Fig. 4. 6^y as a function of tp with = 109.50° and 6 = 6.00° 

Fig. 5. (Color online) Elliptical beams propagation with different spatial widths of the input 
fields: {a)wx : Wy = 1 : 1, (b) : Wy = 1.4 : 1, (c) : Wy = 2 : 1, where Wy = 10.0/im. 
The contour plots of |A(x, y, z)\ are taken at ,^ = 10. The 6 coefficients are, for the mode 1 
beam: 5^ = 1.549 x lO'^, Sy = 6.557 x lO'^, S^^ = -3.906 x lO'^m, 5yy = -3.486 x lO'V 
S^y = -2.095 X lO-^m; for the mode 2 beam: 4 = -1.585 x 10^^^ 6y = 3.745 x lO^^^ 
5xx = -3.565 X 10~^m, 6yy = -3.855 x 10"^m, S^y = 1.427 x 10~^m, respectively. 

Fig. 6. (Color online) Similar to figure 5 except with ip = 69.48°. The 6 coefficients are, 
for the mode 1 beam: 6^ = 1.682 x 10"^ 6y = 1.349 x 10"^ = -3.992 x 10~^m, 
6yy = -3.400 X 10-^m; for the mode 2 beam: 6^ = -4.639 x 10-^ 6y = 4.067 x 10^^ 
Sxx = -3.507 X 10-®m, 6yy = -3.913 x lO'^m, respectively. 

Fig. 7. (Color online) Contour of magnitude of the slowly varying amplitude |A(r)|, for 
an initial Gaussian beam with A = 514nm,?/' = 109.50°, 6 = 6.00°, and ip = 47.00°. Spatial 
widths Wx = 12.5/im and Wy = lOyum in the frame moving with the center of the beam for 

(a)e = 0, (b)e = 0.734, (c)e = 1.277, (d)e = 1.915. 

Fig. 8. (a)Rotation angle <y9 as a function of ^. (b) major to semi major axes ratio 77. Solid 
curve for the mode 1 beam and dashed curve for the mode 2 beam. 

Fig. 9. (Color online) Propagation of mode 1 beam in the assumed metamaterial. Spatial 
widths are Wx = 5/im and Wy = Ifim in the frame moving with the center of the beam for 
(a)e = 0, (b)e = 0.160, (c)e = 0.315, (d)^ = 0.478. 
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Fig. 2. The " x- convent ion" of Euler angles. The first rotation is by angle (f) 
about the 2;'- axis, the second rotation is by an angle ^ G [0, vr] about the former 
X-axis (now oa), and the third rotation is by an angle ip about the 2;-axis. a, 
(3 and 7 are the direction angles of ^-axis in the principle coordinate system 
{x',y',z'). 
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Fig. 3. Different values of ^/^ correspond to different input fields for the ex- 
traordinary beam in uniaxial crystal. For = 7r/4, the input field Eg. (39) 
corresponds to the dashed ellipse and this beam will rotate during propaga- 
tion; for ■?/; = 7r/2, the input field Eg. (39) corresponds to the solid ellipse and 
this beam will not rotate during propagation 
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Fig. 4. b^y as a function of i) with = 109.50° and Q = 6.00° 
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Fig. 5. (Color online) Elliptical beams propagation with different spatial widths 
of the input fields: (a)^^; : Wy = 1 : 1, (b) Wy = 1.4 : 1, (c) '■ Wy = 2 : 1, 
where Wy = 10.0/im. The contour plots of \ A{x, y, z) \ are taken at = 10. The 
6 coefficients are, for the mode 1 beam: 6x = 1.549 x 10~^, 6y = 6.557 x 10~^, 
= -3.906 X 10-*m, 6yy = -3.486 x lO'^m, S^y = -2.095 x lO^^m; for the 
mode 2 beam: 6^ = -1.585 x 10-^ 5y = 3.745 x lO'^, = -3.565 x lO'V 
6yy = -3.855 X 10"^m, d^y = 1.427 x 10"^m, respectively. 
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Fig. 6. (Color online) Similar to figure 5 except with tp = 69.48°. The 6 co- 
efficients are, for the mode 1 beam: 6^ = 1.682 x 10"^, 6y = 1.349 x 10""^, 
6^^ = -3.992 X 10-^m, 6yy = -3.400 x lO'^m; for the mode 2 beam: 4 = 
-4.639 X 10-^ (5y = 4.067x10-3, = -3.507 x 10" V <^?;y = -3.913 xlO'V 
respectively. 




Fig. 7. (Color online) Contour of magnitude of the slowly varying amplitude 
|A(r)|, for an initial Gaussian beam with A = 514nm,?/; = 109.50°, 6 = 6.00°, 
and ip = 47.00°. Spatial widths Wx = 12.5/im and Wy = 10/im in the frame 
moving with the center of the beam for (a).^ = 0, (b).^ = 0.734, (c).^ = 1.277, 
(d)e = 1.915. 
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Fig. 8. (a)Rotation angle ip as a function of (b) major to semi major axes 
ratio rj. Solid curve for the mode 1 beam and dashed curve for the mode 2 
beam. 




Fig. 9. (Color online) Propagation of mode 1 beam in the assumed metama- 
terial. Spatial widths are Wx = 5/im and Wy = 1/im in the frame moving with 
the center of the beam for (a)^ = 0, (b)^ = 0.160, (c)^ = 0.315, (d)^ = 0.478. 
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